##########################################################################################

library(data.table)
library(optparse)
library(ArchR)
library(ggthemes)
library(dplyr)
library(ggsci)
library(ggplotify)
library(parallel)

##########################################################################################
option_list <- list(
    make_option(c("--scriptPath"), type = "character"),
    make_option(c("--comine_data_all_file"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){

    ## 所有的细胞的,计算maxdelt
    comine_data_all_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata"

    ## 既往研究整理的代码
    scriptPath <- "~/20231121_singleMuti/scripts/scScalpChromatin"

    ## 认为关键的TF
    gene <- "DMRT1"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/report_tf"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

scriptPath <- opt$scriptPath
comine_data_all_file <- opt$comine_data_all_file
gene <- opt$gene
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

grp_order2 = c("SSC",
"Differenting&Differented SPG",
"Leptotene",
"Zygotene",
"Patchytene",
"Diplotene",
"Early stage of spermatids",
"Round&ElongateS.tids",
"Sperm",
"Leydig cells",
"Myoid cells",
"Pericytes",
"Sertoli cells",
"Endothelial cells",
"NKT cells",
"Macrophages")

use_colors <- c(pal_npg("nrc")(10) , pal_jco("default")(6))
#names(use_colors) <- unique(scrnat$cell_type)
names(use_colors) <-c( "Myoid cells","Leydig cells" ,          
"Endothelial cells","Zygotene",         
"Round&ElongateS.tids","Patchytene",
"SSC","Sperm" ,                 
"Diplotene","Early stage of spermatids", 
"Leptotene","Sertoli cells",            
"Macrophages","Differenting&Differented SPG",
"Pericytes","NKT cells")

#names(use_colors) <- paste0( as.numeric(factor( names(use_colors) , levels = grp_order2 , order = T )) , names(use_colors))

###########################################################################################
## 导入数据
b <- load(comine_data_all_file)
## testis_combined_peak_combineRNA
atac_proj <- testis_combined_peak_combineRNA

## 细胞类型排序
#atac_proj@cellColData$cell_type <- factor( atac_proj@cellColData$cell_type , levels = grp_order2 , order = T )
#atac_proj@cellColData$cell_type_plot <- paste0(as.numeric(atac_proj@cellColData$cell_type) , atac_proj@cellColData$cell_type)

###########################################################################################
## 已发表文献写好的脚本
if(1!=1){
    source(paste0(scriptPath, "/plotting_config.R"))
    source(paste0(scriptPath, "/misc_helpers.R"))
    source(paste0(scriptPath, "/matrix_helpers.R"))
    source(paste0(scriptPath, "/archr_helpers.R"))
    source(paste0(scriptPath, "/GO_wrappers.R"))
}

###########################################################################################

corrCutoff <- 0.5
label_genes <- gene

p2g <- getPeak2GeneLinks(
  ArchRProj = atac_proj,
  corCutOff = corrCutoff,
  resolution = 100,
  returnLoops = TRUE
)

p <- plotBrowserTrack(
    ArchRProj = atac_proj, 
    geneSymbol = label_genes, 
    groupBy = "cell_type", 
    useGroups = grp_order2,
    pal = use_colors,
    upstream = 250000,
    downstream = 250000,
    loops = p2g
)

out_file <- paste0( out_path , "/" , gene , "_plotBrowserTrack.pdf" )
pdf(out_file,width = 6, height = 7)
#grid::grid.newpage()
#grid::grid.draw(p[[gene]])
#p <- as.ggplot(p[[gene]])
print(as.ggplot(p[[gene]]))
dev.off()

if(1!=1){
    ###########################################################################################
    # Tracks of genes:
    # (Define plot region based on bracketing linked peaks)
    promoterGR <- promoters(getGenes(atac_proj))

    # markerGenes <- c("IL21", "RUNX3")
    mPromoterGR <- promoterGR[promoterGR$symbol %in% label_genes]

    # Restrict to only loops linking genes of interest (full project loops)
    plotLoops <- getPeak2GeneLinks(atac_proj, corCutOff=corrCutoff, resolution = 100)[[1]]
    sol <- findOverlaps(resize(plotLoops, width=1, fix="start"), mPromoterGR)
    eol <- findOverlaps(resize(plotLoops, width=1, fix="end"), mPromoterGR)
    plotLoops <- c(plotLoops[from(sol)], plotLoops[from(eol)])
    plotLoops$symbol <- c(mPromoterGR[to(sol)], mPromoterGR[to(eol)])$symbol
    plotLoops <- plotLoops[width(plotLoops) > 100]

    # Bracket plot regions around loops
    plotRegions <- lapply(label_genes, function(x){
      gr <- range(plotLoops[plotLoops$symbol == x])
      lims <- grLims(gr)
      gr <- GRanges(
          seqnames = seqnames(gr)[1],
          ranges = IRanges(start=lims[1], end=lims[2])
        )
      gr
      }) %>% as(., "GRangesList") %>% unlist()

    # Create copy of individual project plot loops
    sub_plot_loop_list <- list()
    for(pn in names(plot_loop_list)){
      subPlotLoops <- plot_loop_list[[pn]]
      sol <- findOverlaps(resize(subPlotLoops, width=1, fix="start"), mPromoterGR)
      eol <- findOverlaps(resize(subPlotLoops, width=1, fix="end"), mPromoterGR)
      subPlotLoops <- c(subPlotLoops[from(sol)], subPlotLoops[from(eol)])
      subPlotLoops$symbol <- c(mPromoterGR[to(sol)], mPromoterGR[to(eol)])$symbol
      sub_plot_loop_list[[pn]] <- subPlotLoops[width(subPlotLoops) > 100]
    }

    p <- plotBrowserTrack(
        ArchRProj = atac_proj, 
        groupBy = "cell_type", 
        useGroups = grp_order2,
        pal = use_colors,
        plotSummary = c("bulkTrack","featureTrack","loopTrack","geneTrack"), # Doesn't change order...
        sizes = c(7, 0.2, 1.25, 2.5),
        geneSymbol = gene, 
        #region = plotRegions, 
        upstream = 250000,
        downstream = 250000,
        tileSize=500,
        minCells=200,
        loops = p2g
    )

    out_file <- paste0( out_path , "/" , gene , "_plotBrowserTrack.pdf" )
    pdf(out_file,width = 6, height = 7)
    #grid::grid.newpage()
    #grid::grid.draw(p[[gene]])
    #p <- as.ggplot(p[[gene]])
    as.ggplot(p[[gene]])
    dev.off()
}